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We present a modified distance measure for use with distance transforms of anti-aliased, area sampled 
grayscale images of arbitrary binary contours. The modified measure can be used in any vector-propaga- 
tion Euclidean distance transform. Our test implementation in the traditional SSED8 algorithm shows a 
considerable improvement in accuracy and homogeneity of the distance field compared to a traditional 
binary image transform. At the expense of a 10x slowdown for a particular image resolution, we achieve 
an accuracy comparable to a binary transform on a supersampled image with 16x16 higher resolution, 
which would require 256 times more computations and memory. 

© 2010 Elsevier B.V. All rights reserved. 



1. Introduction 

Traditionally, a distance transform (DT) is defined on a binary 
image as the mapping from each foreground picture element (pix- 
el) to its distance from the closest background pixel. Distance func- 
tions and the corresponding DTs are important tools in image 
processing and analysis. An early paper on the subject is Rosenfeld 
and Pfaltz (1966), where a two-scan sequential algorithm was pre- 
sented to compute the DT based on the city-block and chessboard 
distance functions. These distance functions are path-generated, 
since the distance between two pixels is defined as a shortest path 
between the pixels. Several extensions of path-generated distance 
functions have been suggested, see, e.g., Borgefors (1996), 
Rosenfeld and Pfaltz (1968) and Strand (2007). 

In Danielsson (1980), an algorithm based on vector propagation 
for computing the Euclidean DT is presented. The basic idea is to 
assign the vector from the closest background pixel to each fore- 
ground pixel. This approach is also considered in, e.g., Mullikin 
(1992) and Ragnemalm (1993). 

There are many other ways to compute (approximate) Euclid- 
ean DTs. The most common alternative methods are: 

• The fast-marching method, (Sethian, 1999). This frequently used 
method, which can be formulated as a level-set method, 
approximates a differential equation (the Eikonal equation) to 
model an evolving wavefront of constant velocity. Scalar values 
are propagated to compute the arrival time (Euclidean distance) 
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of the wavefront. In Kimmel et al. (1996), sub-pixel distance 
values are computed by interpolation of distance values from 
distance values at adjacent grid points. 

• Methods based on dimensionality reduction, (Felzenszwalb and 
Huttenlocher, 2004; Maurer et al., 2003; Saito and Toriwaki, 
1994): Since the squared Euclidean distance value is separable, 
the DT can be computed by one-dimensional scans along each 
principal direction. 

• Grayscale morphology: The (squared) Euclidean DT can be 
obtained by dilating the binary image with a (large) structuring 
element. Efficient algorithms are obtained by decomposing the 
structuring element, see Huang and Mitchell (1994), Mehnert 
and Jackway (1999) and Shih and Mitchell (1992). 

2. Overcoming the constraints of a binary image 

The traditional restriction of distance transforms to binary 
images, where each pixel belongs to either the foreground or the 
background, results in a well defined discrete problem, and in 
many applications the constraints of a binary input image can be 
tolerated. However, a binary image is spatially discretized to a reg- 
ular grid, and therefore unable to represent general binary shapes 
correctly. Edges which are not aligned with the sampling grid are 
discretized with strong stair-step edge aliasing artifacts, referred 
to in the computer graphics literature as as "jaggies". A binary 
two-dimensional shape with crisp edges of arbitrary position and 
orientation contains significant high spatial frequency components 
and is not band limited. It is not well represented by point sam- 
pled, binary-valued pixels and can not be well reconstructed from 
such data (Kothe, 2008; Schilling, 2001). 

A binary image conveys no information on the sub-pixel posi- 
tion of the edge of the sampled shape, as each sample is considered 
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to be either fully inside or fully outside the contour. Therefore, dis- 
tance vectors for a binary DT are constrained to an integer grid. Fol- 
lowing Sladoje and Lindblad (2009) and common practice in 
computer graphics, we use a model where the image is instead 
area sampled, whereby the intensity value in a pixel represents 
the area coverage, the fraction of the pixel rectangle that is inter- 
sected by an object. The result is that pixels intersecting the edge 
of a binary contour attain grayscale values, whereby edge aliasing 
is greatly reduced. The limited sub-pixel information made avail- 
able through the area coverage can be used to improve precision 
in, e.g., estimators for edge length and position (Verbeek and van 
Vliet, 1993), arc length and area (Eberly and Lancaster, 1991), 
and reconstruction of the digitized object (Kiryati and Bruckstein, 
1988). Using the area sampling model, Sladoje and Lindblad 
(2009) gives a method that provides error-free measurements of 
the length of straight boundary segments in the case of non-quan- 
tized pixel values. 

Real world digital imaging sensors consist of finite-area sensor 
elements, and the sampling performed by such devices is inherently 
an area sampling, with light being integrated over the sensor ele- 
ment convolved with the point spread function of the optics (Chen 
et al., 2009). While the actual sampling kernel of a real imaging de- 
vice does not have a crisp edge or a rectangular shape, it is certainly 
much closer to an area sampling than to a point sampling. In syn- 
thetic, rendered images, point sampling causes aliasing of high fre- 
quency features in the image, e.g. edges. Multi-sampling or 
analytic anti-aliasing strategies are employed to effectively replace 
the point sampling with an approximation of area sampling, bring- 
ing the sampling closer to what a physical imaging device would do. 

In the novel DT method presented in this paper, the distance 
vectors point not to the center of the closest pixel in the spatially 
discretized representation of the object, but instead to the closest 
position on the estimated border of the spatially continuous object. 
This extension of the traditional definition of the distance in DT 
algorithms makes the distance vectors unconstrained by the dis- 
crete sample grid. The position of the border is estimated by using 
the area coverage of a single pixel straddling the edge. An im- 
proved estimate additionally uses a small region of adjacent pixels 
for better accuracy near the edge. Scalar propagation methods 
(listed above) can often be extended to use area sampled grayscale 
images instead of point sampled binary images as input. Grayscale 
valued input images are considered for the dimensionality reduc- 
tion method in Felzenszwalb and Huttenlocher (2004). Another 
scalar propagation method with sub-pixel accuracy is presented 
in Krissian and Westin (2005), where interpolation of the com- 
puted distance values in a level set method is used to get a better 
approximation of the wavefront that corresponds to the evolving 
level set. In Krissian and Westin (2005), the level set is approxi- 
mated by interpolation using distance values at grid points in a 
narrow strip and, similarly to our method, the sub-pixel accuracy 
is obtained by using the distance to the underlying continuous 
curve instead of the distance to the discrete grid points. 

Under the assumption that the edge is locally straight, the pixel 
coverage and the edge direction together determine the distance 
from the pixel center to the edge. In this paper, the vector to the 
closest pixel intersecting the edge is adjusted by a small sub-pixel 
distance computed from the approximate direction of the edge and 
the area coverage of the edge pixel. The result is a sub-pixel accu- 
racy extension suitable for use with any vector propagation meth- 
od for DT computation. 



3. Applications 

There are a number of image processing applications where the 
proposed method can be used to provide better results than a binary 
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DT. By modifying the definition of distance to measure the distance 
to the ideal boundary between foreground and background, signed 
distance fields created by combining two DT passes over the image, 
one on the image and one on its inverse, will not exhibit a problem- 
atic discontinuity at the edge, and therefore morphological opera- 
tions like erosion and dilation with small distances of only a few 
pixels will be significantly more accurate than by using a binary 
DT. Euclidean erosion and dilation can be performed with arbitrary 
fractional distances, even less than one pixel. This is not possible 
with a binary transform (Ragnemalm, 1 992). Such sub-pixel erosion 
and dilation of anti-aliased images are useful for e.g. digital photoli- 
thography, where there is a need to compensate for imperfections in 
the exposure and etching process which yield an erosion or dilation 
of the final outlines (Wei and Brainard, 2009). Level sets created from 
DT images with our method will be more accurate, in particular near 
the edge, even from lower resolution input data. Accurate distance 
maps for synthetic images can be created from regular anti-aliased 
images instead of from special, high resolution binary images, which 
are cumbersome to create and difficult to process (Green, 2007). 
Accurate distance maps can also be created from (approximately) 
area sampled digitized images of real world objects, better utilizing 
the information in images where the spatial resolution is limited by 
physical constraints. 

4. Anti-aliased distance function 

A binary distance transform is defined for images with only two 
pixel values: each pixel is considered to belong either entirely to 
the foreground or entirely to the background. To lift this restric- 
tion, we make the following definitions: Assume a piecewise con- 
stant binary function B(x,y)e[0, 1} defined over a continuous 
plane (x,y). B{x,y) contains a set of compact regions of arbitrary 
shape and with discontinuous step edges. By arbitrary choice, 
B = 1 represents foreground and B = 0 the background. (The oppo- 
site convention is also quite common.) We create a square grid 
of pixels and assign a value a, 0 C a ^ 1, to each pixel based on 
the fraction of the pixel that is intersected by the foreground, as 
illustrated in Fig. 1(a) and (b). The grayscale values in the input im- 
age are used to approximate the location of the object boundary. 
The difference between the traditional binary distance measure 
and our modified measure is illustrated in Fig. 1(c). 

Our objective is to find the distance to the closest ideal edge 
from a discrete set of area sampled pixels by analyzing only the 
pixel values. If we know the orientation of the underlying edge 
and assume that its local radius of curvature is significantly greater 
than the pixel size, we can compute the sub-pixel edge position 
within the pixel from the area coverage a. Denote the (normalized) 
edge direction cp, which defines g = (g x ,g y ) = (coscp,sincp). The edge 
direction cp determines the areas a A and a 2 , the distance values d] 
and d 2 illustrated in Fig. 2. 

The signed distance from the edge to the center of the pixel is de- 
noted df, see Fig. 2, right. Under the assumption that the edge of the 
object covering the pixel is straight, there is a relation between df, 
the area coverage a, and the edge direction cp. For all orientations 
of the edge, a = 0.5 means df = 0, and for g x = 0 or g y = 0, the relation 
between df and a is a simple linear ramp: 

df = 0.5 -a. (1) 

Note that df = 0 represents an edge that passes through the center of 
the pixel. For edges in other than axis-aligned directions, this 
approximation can be off by as much as ±(\/2 - l)/2 but is still use- 
ful as a simple, first order approximation. If the edge direction is 
known, df can be more accurately expressed as a function of {g x ,g y ) 
and a, see Fig. 2. Using elementary trigonometry, we can derive the 
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Fig. 1. (a): Pixel representation of a disc by area coverage. True contour shown by dashed line, (b): Detail view of upper right corner of (a), (c): The difference between the 
traditional integer distance definition (dashed arrow) and the definition proposed in this paper (solid arrow). 
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Fig. 2. Notation for working out a as a function of df and vice versa. See Eq. 2 for details. 



areas di and a 2 of the left-hand side of Fig. 2 and the distances di, d 2 
in the middle: 

a, = tan cp = g y /g x , 
a 2 = 1 — 20i, 
d, = sin<p = g y , 

d 2 = (l/\/2)sin(7r/4- cp) = (1/V2) sin(7t/4 - arcsing y ). (2) 

The continuous and monotonous function a{d f ) for a given edge 
direction ip can then be written as: 



0. 

(d,+d 2 + d f ) 2 (a,/d]), 
a= \a^ + (a 2 /2)(l +d f /d 2 ), 

1 -(d, +d 2 -d / ) 2 (a,/d'), 
1, 



-d] - d 2 , 
-d] - d 2 d ; < - 
-d 2 < d < d 2 , 
d 2 < d f < di +d 2 , 
df < di +d 2 . 



(3) 



In Eq. 3, we can see that, for an arbitrary edge direction <p, the 
area coverage a for a pixel as a function of d f starts out as a parab- 
ola (proportional to d 2 ) when the edge is in the corner area a A 
shown to the left in Fig. 2, followed by a linear ramp (proportional 
to df) within the center area a 2 and a final parabola in the opposite 
corner. The function is antisymmetric around d/=0 and a = 0.5, 
such that a(-df) = 1 - a{d f ), and symmetric in tp such that it is en- 
ough to look at angles in the first octant g x 3= g y 3= 0. Other cases 
can be evaluated by changing the sign of g x and/or g y and if neces- 
sary swapping g x and g y . 

For our application, we want to compute df as a function of a. 
For a fixed edge direction ip, this is the inverse of Eq. 3. The inverse 
df[a) is defined for 0 C a ^ 1: 



(g x +g y )/2- ^2g x g y a. 
(0.5 - a)g x , 

~(g x +S y )/2 + J2g x g y {\ - a) 



0 < a < a,. 

a\ < a < 1 — ai, 

1 - a! < a < 1. 



(4) 



A plot of this function for three different edge directions is pre- 
sented in Fig. 3. The plot also shows that the linear function 
df= 0.5 - a, which is valid for ip = 0 and ip = 7i/2, is a reasonable 
approximation for any edge direction cp. 



0.8 t d, 




Fig. 3. Plot of d/a) for cp = 0°, 20°, 45°. 



5. The anti-aliased vector propagation method 

In vector propagation distance transform algorithms, the vec- 
tor (d x ,dy) to a candidate for the closest edge is known, and a 
distance function is computed from the vector. A binary Euclid- 
ean distance transform uses the distance function d t = W 1 d 2 x + d y . 
A simple but reasonable assumption is that the edge of the ob- 
ject is approximately orthogonal to the vector (d x ,d y ). This 
assumption does not hold true for small d, where only a few dis- 
crete directions are possible, but it is increasingly accurate for 
pixels further away from the edge, and asymptotically exact 
for distant pixels. 

In the anti-aliased vector propagation method, d,- is adjusted 
using d^ from Eq. 4 above: 



d = di + d f 



(5) 
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A vector propagation distance transform using the modified dis- 
tance measure d requires more computations than merely comput- 
ing d, for each candidate tested, but because df is bounded, an 
early-out test on d, could avoid many computations of d f in a 
sequential implementation. Note also that df is both bounded and 
smooth and can be tabulated with good accuracy using a small 
2D table and bilinear interpolation if speed is important. 

5. J. Further improvements 

At edge pixels, the distance vector points to the same pixel and 
does not contain any information on the direction of the edge. Near 
the edge, the integer vector from a pixel to the closest edge pixel 
does not align well with the true direction of the edge. For greater 
accuracy at those pixels, the local direction of the edge was esti- 
mated by two 3x3 gradient filter kernels around each edge pixel. 
At the edge, the estimated local gradient was always used for 
(gx,g y ) in Eq. 4. Near the edge, a distance measure according to 
Fig. 4 was computed based on the local gradient at the edge and 
used instead of d, + df, but only if it resulted in a well defined 
orthogonal distance to a linear edge with the same gradient, as ex- 
plained by Eq. 6. 

In our experiments, using the local gradient always improved 
the accuracy for edge pixels, where no other edge direction infor- 
mation was available. Away from the edge, using the local gradient 
sometimes increased the error but still improved the average re- 
sult within one or two pixels distance from the edge. Beyond that, 
the gradient information was frequently less accurate than the 
integer distance vector. Therefore, the local gradients were used 
to improve the result only at the edge and very near it, where most 
of the larger absolute errors were located. 

(g x ,g y ) is local gradient at edge pixel 
(d x ,d y ) is vector to the edge pixel 
compute d; from (4) using (g x ,g y ) 
t = dj sin p = g x d y - g y d x (cross product, |g| = 1 ) 

« = ~ d fS x + tg y 
v = -d f g y - tg x 

if\u\ < 0.5 and \v\ < 0.5 ("hit inside pixel") 



d= ^{dx + uf + (dy+ vf 
else 

compute dj from (4) using (d x , d y ) 



d=Jdt + dl + d, 



(6) 



5.2. Results 



We implemented our modified distance measures in a classic 
sequential Euclidean distance transform (EDT) vector propaga- 
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tion algorithm using eight neighbours (8SED), see Danielsson 
(1980). The implementations were not optimized for speed, as 
we were only interested in showing the performance relative 
to a classic binary DT. The exact same propagation code was 
used for both the reference binary DT and the improved anti-ali- 
ased DT, so the performance figures in Table 1 are good indica- 
tors of the relative slowdown due to using the new floating- 
point distance measure in the propagation loop instead of the 
integer-only evaluation of df for the binary DT. The performance 
figures show that the much improved accuracy comes at a rea- 
sonable cost in execution time, even without any early-out tests 
or look-up table speedups. The exact slowdown for the last 
improvement (Eq. 6) depends somewhat on the amount of edges 
in the image, as the extra effort is spent only on pixels near 
edges. However, the main algorithm requires several sweeps 
and updates across the image, and the final adjustment of dis- 
tances using the pre-computed local gradient information consti- 
tuted only a minor part of the computations. In a worst case 
scenario where all pixels in the image were close enough to 
an edge to be post-processed in this manner, the relative slow- 
down was only around 3%. 

For experimental evaluation of the performance of our algo- 
rithm, we used a large set of test images of different complexity. 
For this presentation, we choose to focus on one carefully chosen 
representative contour, shown in Fig. 5. This contour contains 
edges in all directions with a wide range of positive and negative 
curvatures, and constitutes a difficult and representative test 
case. Using our modified distance measure and our proposed 
algorithm resulted in significant improvement in accuracy over 
a traditional binary DT, as shown in Fig. 6. The error is greatly 
reduced even with the simple linear approximation of Eq. 1 
(Fig. 6, lower left), and it is further reduced with the more accu- 
rate measure of Eq. 4 taking the vector direction into consider- 
ation (Fig. 6, lower middle). The remaining errors near edges 
are significantly reduced by using gradient filters to supply an 
improved estimate of the edge direction, as presented in Eq. 6 
(Fig. 6, lower right). 

The average absolute error of our final result, the lower right 
image of Fig. 6, is 0.02 pixels. This is significantly better than the 
upper middle image which was computed by a binary EDT on an 
image of 8 x 8 higher resolution, and comparable to the upper 
right image using a binary EDT of 16 x 16 higher resolution. 99% 
of the pixels in our result are accurate to within ± 0.2 pixels. A 
few isolated pixels exhibit absolute errors of up to 0.5. These worst 
case errors occur where the closest contour has a radius of curva- 
ture smaller than one pixel, which invalidates our assumption of a 
locally straight edge. A region of such worst case errors can be seen 
extending from the pointed left fin of the dolphin countour. This 
single point of strong curvature is in fact the main source of error 
in this test image. 

The anti-aliased EDT yields a smoother gradient field, as illus- 
trated in Fig. 8 by the divergence of the gradient of the distance 
field. While the modulus of the gradient is quite close to 1 at most 
positions even for a binary EDT, the gradient direction deviates sig- 
nificantly from the true direction in an unpredictable manner. 
There are noticeable directional discontinuities in the computed 



Table 1 

Relative performance of traditional vs. new measures. Image size 2048 x 2048 pixels, 
test pattern according to Fig. 5, CPU Intel Core2 Duo 2.4 GHz, single thread execution. 
First row is execution time, second row is relative slowdown. 



Fig. 4. Notation for determining the true orthogonal distance d to an edge, given d x , 
d y , g x , g y , df. See also Eq. 6. If the hit point (u, v) (black circle) is inside the pixel, the 
true distance d is used instead of d, + d f . 
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Fig. 5. (left) Test contour for the experiments, (middle) 256 x 256 pixels area-sampled rendering of the contour, (right) Euclidean distance transform of the middle image. 





Fig. 6. Error images for 256 x 256 pixels EDT results: transform result minus true distance. Zero is medium gray, white is + 1, black is -1. Top row: binary EDT, subsampled 
binary EDT from 8x8 higher resolution, subsampled binary EDT from 16 x 16 higher resolution. Bottom row: Approximate AA EDT according to Eq. 1, AA EDT according 
to Eq. 4, further improved AA EDT according to Eq. 6. Errors are discussed in Fig. 7 and in the text. 
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Fig. 7. Histograms for the respective error images in Fig. 6. Systematic errors are eliminated and accuracy is improved by our algorithms. Our final result (lower right) is 
comparable to using a binary EDT of 16 x 16 higher resolution (upper right). For details, see the text. 



distance field even at large distances, as can be seen in Fig. 8, left. 
These artifacts are reduced to the limit of extinction by using an 
anti-aliased grayscale input image and our modified distance 
transform algorithm, see Fig. 8, right. 



6. Conclusion 

The distance field generated by our improved algorithm has 
good accuracy even near contour edges. This makes it well suited 
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Fig. 8. Divergence of the discrete, filter-estimated gradient of the distance field. White is positive, black is negative, medium gray is zero. Left: binary EDT. Right: our anti- 
aliased EDT. 



to performing erosion and dilation of contours with fractional and 
sub-pixel distances. 

A distance field with sub-pixel accuracy is also well suited for 
detailed reconstruction of sharp boundaries of arbitrary shape 
from a sampled image, while a binary transform is inadequate 
for this purpose. This reconstruction problem, in fact, was our 
motivation for developing this transform. 

Furthermore, the improved smoothness of the gradient of the 
generated distance field, as demonstrated by Fig. 8, makes the dis- 
tance field much better suited for generating level sets, height 
fields, bump maps, bevels and other data structures for computer 
graphics where the visual result is dependent on differentials of 
the distance field. 

In theory, the anti-aliased distance transform presented here 
can be replaced by a traditional binary transform on a higher res- 
olution image, followed by subsampling. This is what we used for 
the "ground truth" in our experiments to evaluate the accuracy. 
However, even when a higher resolution binary image can be ac- 
quired (which is typically not the case for real world input data), 
the return in edge position accuracy for each unit of input data is 
much greater with antialiased grayscale images than with binary 
images. A linear edge rendered with ideal area sampling and quan- 
tized to 8 bits per pixel can theoretically convey the edge position 
with an accuracy of around 1/255 of a pixel, somewhat depending 
on the edge position and orientation. If a similar minimum accu- 
racy were to be computed from binary input, the binary image 
would need about 255 x 255 as many 1-bit pixels, or 255 x 255/ 
8 rs 8000 times as much data. Even though our results fall short 
of this theoretical maximum, we show an accuracy equivalent to 
a 16 x 16 higher resolution binary EDT, which translates to around 
30 x savings in execution time and memory requirements. We con- 
clude that it is clearly worthwhile to consider the extra accuracy 
which can be achieved in a distance transform by utilizing gray- 
scale information where available. 

7. Source code 

For easy evaluation and testing, we have wrapped our anti-ali- 
ased Euclidean distance transform implementations in MEX files 
for Matlab or GNU Octave. The source code is published under 
the GPL and can be downloaded from: http://www.itn.liu.se/ 
~stegu/edtaa/. 
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